Germline genetic regulation of the colorectal tumor immune microenvironment

Objective To evaluate the contribution of germline genetics to regulating the briskness and diversity of T cell responses in CRC, we conducted a genome-wide association study to examine the associations between germline genetic variation and quantitative measures of T cell landscapes in 2,876 colorectal tumors from participants in the Molecular Epidemiology of Colorectal Cancer Study (MECC). Methods Germline DNA samples were genotyped and imputed using genome-wide arrays. Tumor DNA samples were extracted from paraffin blocks, and T cell receptor clonality and abundance were quantified by immunoSEQ (Adaptive Biotechnologies, Seattle, WA). Tumor infiltrating lymphocytes per high powered field (TILs/hpf) were scored by a gastrointestinal pathologist. Regression models were used to evaluate the associations between each variant and the three T-cell features, adjusting for sex, age, genotyping platform, and global ancestry. Three independent datasets were used for replication. Results We identified a SNP (rs4918567) near RBM20 associated with clonality at a genome-wide significant threshold of 5 × 10− 8, with a consistent direction of association in both discovery and replication datasets. Expression quantitative trait (eQTL) analyses and in silico functional annotation for these loci provided insights into potential functional roles, including a statistically significant eQTL between the T allele at rs4918567 and higher expression of ADRA2A (P = 0.012) in healthy colon mucosa. Conclusions Our study suggests that germline genetic variation is associated with the quantity and diversity of adaptive immune responses in CRC. Further studies are warranted to replicate these findings in additional samples and to investigate functional genomic mechanisms. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10295-1.


Background
In colorectal cancer (CRC), a strong tumor infiltrating lymphocyte (TIL) response is an established prognostic indicator for better disease-specific survival and a predictive factor for response to checkpoint inhibitor immunotherapy [1][2][3][4][5].However, the extent and diversity of T cell responses in the CRC tumor microenvironment are highly variable with mostly uncertain drivers.Beyond microsatellite instability (MSI) and tumor mutational burden, the factors contributing to heterogeneity in adaptive immune responses, including abundance and clonality of T cells, across tumors remain largely unexplored.Germline variation, particularly in genes important for immune cell differentiation and function, may contribute to the strength and quality of adaptive immune responses mounted against colorectal tumor cells.In turn, this variability may influence both the risk of developing cancer and subsequent clinical outcomes following a CRC diagnosis via altered immunosurveillance.
Indeed, examples of germline genetic variation modifying the diversity of constitutional immune environments and host immune responses to cancer are numerous.Single nucleotide polymorphisms (SNPs) in the Major Histocompatibility Complex locus containing HLA genes are highly associated with risk of Hodgkin lymphoma [6].Further, pro-inflammatory cytokines have been strongly associated with African germline genetic ancestry proportion among African American women attributable to a Duffy-null allele [7].Recently, three separate pancancer investigations utilizing The Cancer Genome Atlas (TCGA) data revealed multiple inherited genetic factors to be associated with tumor-immune compositional and functional features of cancer, including CRC [8][9][10].However, each included only a limited number of CRCs and relied upon an estimation of immune infiltration from transcriptomic data.In other investigations specific to CRC, common variants in immune-related genes and pathways have been linked to CRC risk and outcomes [11][12][13][14].For example, common SNPs in the T regulatory cell pathway gene, TGFBR3, have been associated with CRC survival [15].Finally, evidence of strong genetic susceptibility variants for autoimmune and other immunerelated diseases strengthens the premise for germline regulation of adaptive immune responses.To build upon prior studies, we conducted the largest and most comprehensive study to date of common germline genetic variability in relation to features of T cell responses in CRCs from population-based studies.

GWAS of T cell receptor clonality
The discovery GWAS identified 18 SNPs with p < 5 × 10 − 6 residing in 17 genomic regions associated with immuno-SEQ-derived clonality.Of these, six variants at 2q24.1, 4q21.21,6q21, 9p13.3,10q25.2, and 18q21.32remained below this statistical significance threshold upon joint discovery-replication meta-analysis (Table 1; Figure S2,  S3).Although the associations from the joint analyses were mainly driven by the discovery data, the association results from the smaller replication sets still demonstrated the same direction and similar magnitude.The strongest association signal in the discovery phase arose from a common variant located at 10q25.2 [rs4918567; odds ratio (OR): 0.76 (95% confidence interval (CI): 0.69-0.84);p = 3.73 × 10 − 8 ].This variant remained associated with clonality in the joint meta-analyses at the genomewide significant threshold with no evidence of heterogeneity (Table 1).eQTL analysis from healthy colon mucosa using BarcUVa-Seq showed that homozygous T/T individuals for rs4918567 had statistically significantly higher expression in ADRA2A (p = 0.012; Figure S4).
A common variant residing in the intronic region of FRAS1 at 4q21.21 [rs4443313; OR: 0.83 (95% CI: 0.78-0.89);p = 5.26 × 10 − 8 ], approaching genome-wide statistical significance, was identified in the discovery phase.However, the strength of the association was weakened in the joint meta-analysis due to a different direction of association in one of the replication sets [P for heterogeneity = 0.055, Table 1].This variant was associated with the differential expression of LINC01094 in eQTL analysis using normal colon mucosa samples (p = 0.0025, Figure S5) and CCNG2 in tumor tissue samples (p = 0.024, Figure S6).
The other association signals at 2q24.1, 9p13.3 and 18q21.32involved common variants.Of note, SNP rs7874748 at 9p13.3 was an eQTL for several genes including SPAG8 (using GTEx/v7 transverse colon), RECK and NPR2 (healthy mucosa from BarcUVa-Seq; Figure S7).rs7874748 was also an eQTL for several genes including RUSC2, CCDC107, CA9, TPM2, TLN1, MSMP, NPR2, OR13J1, and HRCT1 in tumor tissues (using Conclusions Our study suggests that germline genetic variation is associated with the quantity and diversity of adaptive immune responses in CRC.Further studies are warranted to replicate these findings in additional samples and to investigate functional genomic mechanisms.

GWAS of T cell receptor abundance
In the discovery phase, we identified 19 SNPs across 17 genomic regions associated with abundance with p < 5 × 10 − 6 .Of those, 11 SNPs remained with p < 5 × 10 − 6 in the joint meta-analysis (Table 2; Figure S10, S11).The strongest association signal identified in the discovery phase was for a common variant (rs1518405) at 2q31.3 with a p value of 8.67 × 10 − 7 (OR: 1.59, 95% CI: 1.32-1.91).However, no significant results for chromatin interactions or eQTL analysis were found.Of note, rs56148061 at 11q12 was associated with differential expression of SYT7 in eQTL analysis using colorectal tumor samples (p = 0.007, Figure S12).

Subset analysis in MSS tumors
When restricting to MSS participants only, the directions of associations and significance levels remained consistent for the top 19 variants from the GWAS.Five additional SNPs were associated with p < 5 × 10 − 8 specific to MSS tumors were identified (Table S2-S4).

T cell subtype infiltration and top variants from GWAS
Of the 19 top variants from our GWAS analyses, one SNP was significantly associated with cytotoxic T cell Violin plots for best-call genotypes of each SNP and abundance are included in Supplementary Figure S19 infiltration as measured by mIF.rs1518405, with the strongest association for abundance (p = 4.08 × 10 − 7 ), was significantly associated with the infiltration of CD3 + /CD8 + T cells.Individuals who carry an A allele at this locus were associated with higher density of cytotoxic T cells (per allele OR: 1.54, 95% CI: 1.02-2.29,p = 0.0386), showing the same direction of association as we observed for abundance (OR: 1.59; 95% CI: 1.33-1.90).Another variant, rs577783 at 6q23.3, with a significant association in the joint meta-analysis for abundance (p = 1.23 × 10 − 6 ), was associated with the infiltration of memory T cells (OR: 1.10, 95% CI: 1.01-1.20,p = 0.0316) and regulatory T cells (OR: 0.69, 95% CI: 0.49-0.97,p = 0.0351).

Associations between previously reported variants with T cell features
Previously reported germline variants associated with immune infiltration in tumors were examined in the Discovery GWAS and the joint meta-analyses (Table S5).In a large study of the association between variants in immune related disorders and CRC risk, the authors reported two variants, rs11676348 and rs102275 [16].Of those, we detected an association between rs102275, a susceptibility locus for Crohn's disease, and abundance in the Discovery GWAS and joint meta-analyses (p = 0.0205 and 0.0061, respectively; Table S5).Another study demonstrated associations between two germline SNPs, rs3366 and rs4819959, and the amount of follicular helper T cells in all solid tumors using TCGA data [8].These SNPs were not associated with any T cell features in our study.Lastly, of the 1,587 germline variants associated with 33 immune traits [9], we were able to examine 1,323 SNPs in our study but did not observe any associations between these variants and the T cell features in our data (Tables S6).
Common variants associated with CRC risk at a genome-wide significant threshold from previous GWAS were also investigated in relation to the three T cell features assessed in this study (Table S7).Among the 155 risk alleles, we identified 7, 9, and 7 SNPs (for clonality, abundance, and TILs, respectively) associated with T cell outcomes in the joint meta-analysis at a nominal level of statistical significance (p < 0.05).

Discussion
Tumor immune microenvironments differ substantially between patients, yet the regulatory control of this diversity is not well understood.Here, we present evidence from a unique study of well-characterized, populationbased, incident CRCs to explore inherited contributions to T cell responses to CRC.We showed that adaptive immune responses are partially determined by the inherited genome and identified specific candidate loci that contribute to both the intensity and diversity of T cell responses, as reflected by quantifiable measures of tumor infiltrating lymphocytes as well as specific clonal responses.
Our strongest association signal in the GWAS for T cell receptor clonality came from rs4918567 at 10q25.2 in the intronic region of RBM20.RBM20, encoding a protein that binds RNA and regulates splicing, is associated with familial cardiomyopathy [17].A recent study indicates a potential role for RBM20 in cardiovascular complications of diabetes by mediating insulin damage in cardiac tissues [18].eQTL analyses for rs4918567 showed an association between homozygote T/T and higher expression in ADRA2A.ADRA2A, located downstream of RBM20, was associated with neutrophil percentage of leukocytes [19] and type 2 diabetes [20,21].Higher expression of ADRA2A was also associated with breast cancer survival and hypothesized to suppress cell proliferation and invasion through the PI3K/AKT/MTOR pathway in cervical cancer [22,23].Another variant rs4443313, located in the intronic region of FRAS1, was strongly associated with clonality.eQTL analyses for this SNP showed differential expression of CCNG2 in colon tumor tissues.
We also identified a rare variant (rs184508436, MAF = 0.009), residing in the intergenic region of ATG5 and PRDM1, to be strongly associated with clonality.ATG5 encodes a protein involved in autophagic vesicle formation, negative regulation of innate antiviral immune response, and apoptosis [24].A recent pan-cancer study using multiple public databases showed that the expression of ATG5 is associated with tumor immune infiltration in most solid tumors [24].PRDM1 is a protein coding gene, which involves in the pathway of regulation of TP53 expression and degradation.SNPs in the region of PRDM1 have been associated with Crohn's disease and inflammatory bowel disease [25][26][27][28].Although no significant eQTL results were identified for this variant, the potential role of ATG5 in tumor metabolism and immune escape warrants further investigation.
The strongest association signal identified for T cell receptor abundance was rs1518405, located at 2q31.3 in the intergenic region of CWC22.Variants in CWC22 have been associated with various traits, including schizophrenia, body mass index, and prostate cancer [29][30][31].rs1518405 was also significantly associated specifically with the infiltration of CD3 + /CD8 + T cells from our subset analysis using mIF.This preliminary finding provides evidence that the association between rs1518405 and abundance may be particularly driven by the infiltration of cytotoxic T cells.
We identified two variants associated with TILs.rs10982853, located in the intergenic region between DELEC1 and PAPPA at 9q33, was significantly associated with several genes in this region from the eQTL and chromatic interaction mapping analyses.This region has been associated with blood protein levels and lymphocyte counts [19,[32][33][34][35]. PAPPA plays an important role in the activation of the IGF pathway.Overexpression of this protein may enhance IGF receptor signaling and promote tumor growth and invasion [36].The second variant associated with TILs, rs215529, was located between a pseudogene RPL19P1 and CPXM1 at 20p13.FASTKD5 and CPXM1, linked by chromatin interaction mapping, have been previously associated with white blood cell counts and blood protein levels [33,35].Interestingly, a recent study demonstrated that an immune-related gene prognostic index constructed using SFRP4, CPXM1, and COL5A was associated with better survival and better responses to immune checkpoint inhibitor therapy for head and neck squamous cell carcinoma [37].CPXM1 was also included in a tumor microenvironment prognostic signature profile comprised of 11 immune checkpoint genes for advanced-stage serous ovarian cancer [38].
Given that the tumor microenvironments in CRC differ by the MSI status, we analyzed the top variants from our study limited to the subset of patients with MSS tumors.When restricting analyses to patients with MSS tumors only, the top variants from our study showed the similar magnitude of effect sizes and p values remained significant (Table S2-S4).Furthermore, we performed analyses by treating MSI status as a confounder.Additional adjustment for MSI status in the discovery dataset did not yield any notable changes in the top findings (data not shown).
Upon examining germline variants identified in previous reports, we detected an association between rs102275 and T cell receptor abundance in our joint meta-analyses (p = 0.0061, Table S5).This variant at 11q12.2 is a risk SNP for Crohn's Disease and inversely associated CRC risk (OR = 0.92, p = 2 × 10 − 5 ). 16Our study showed that rs102275 may be associated with adaptive immune responses, especially the intensity of T cell infiltration.rs102275 was also located in proximity to rs174537 from 155 variants reported by a prior CRC GWAS (Table S7).rs174537 is significantly associated with CRC risk in individuals of East Asian and European descent (p = 9.22 × 10 − 21 and 7.39 × 10 − 5 , respectively) [39].It is also an eQTL for FADS1 and FADS2.Although we did not observe an association for rs174537, rs102275 is in high LD with rs174537 (R 2 = 0.933; D'=0.996) in European populations.Of note, rs102275 has been associated with several traits related to lipid metabolism and blood metabolite measurement from the NHGRI-EBI GWAS Catalog [40].
Our study is the first and largest population-based study compiling clinical, pathological and epidemiologic data with high-quality immunoSEQ results to investigate the association between germline genetic variation and T cell infiltration in CRC.Our findings suggest that germline genetic variation is associated with the quantity and quality of adaptive immune responses in CRC patients.To our knowledge, this unique study utilizes all existing relevant datasets and provides insight into the role of germline genetic variation in relation to the T cell repertoire, despite the replication dataset sample sizes.Nonetheless, we acknowlege that our study has several limitations.First, because the sample sizes of the three replication datasets were relatively small, the metaanalysis results were primarily driven by the discovery study.Likely due to the small sample sizes of the replication cohorts, most of the suggestive loci (p < 5 × 10 − 6 ) from the discovery data were not independently replicated.Second, we recognize that using the conventional GWAS significance threshold (p < 5 × 10 − 8 ) is a less stringent criterion for a larger number of imputed SNPs, including those with low-to-intermediate frequency, tested in our analyses.Importantly, here we present the combined results from all currently available epidemiological studies with genomic data and T cell immune features in CRC with the goal to highlight suggestive loci for the association between germline variations and T cell features in the CRC tumor microenvironment.Further studies with larger samples are warranted to replicate our findings and to identify the underlying functional genomic mechanisms.Because immune function is a key determinant of immunotherapy outcomes, there are important clinical implications of understanding the factors that influence the robustness and diversity of T cell responses in the CRC tumor microenvironment [41][42][43][44].In CRC, several immune checkpoint inhibitors are FDAapproved for the treatment of metastatic microsatellite instable CRC, a molecular subtype of CRC in which TILs are enriched [43].Our study underscores the importance of future work to determine the role of germline genetic regulation in response to immune checkpoint inhibitors.Germline genetic variation is likely to contribute to the wide range of responses to immune checkpoint inhibitors, and the results presented here provide direction for exploring candidate genes that are likely to regulate the adaptive immune responses in CRC.This presents the possibility that pharmacogenetic variation in the genes

Discovery Phase Study Population
The Molecular Epidemiology of Colorectal Cancer Study (MECC) is a population-based study of incident CRC cases and healthy controls recruited from northern Israel from 1998 to 2017 [45].Cases included patients with pathologically confirmed invasive colorectal adenocarcinomas.Controls are participants without a prior history of CRC selected from the same source population as cases and with individual matching on age, sex, Jewish ethnicity, and primary clinic site.Baseline demographic and clinical characteristics of the CRC cases contributing to this study are described in Table 4.Of note, the distributions of key demographic and clinical characteristics of all MECC cases and the subsets of cases whose tumors underwent pathological review and/or immunosequencing did not substantially differ (Table S1).

Genotyping, quality control, and imputation and PCA analysis for discovery phase (MECC)
Germline DNA was extracted from peripheral blood samples.DNA samples from 10,364 MECC participants (5,581 CRC cases and 4,783 controls) were genotyped using four platforms.All genotype data were cleaned by platform using common quality control metrics at the individual and SNP levels as described previously [46][47][48].485 cases and 498 controls were genotyped in two batches using Illumina HumanOmni2.5 chips, measuring approximately 2.3 million SNPs [47].Batch 1 (384 cases and 143 controls) was genotyped at the Case Western Reserve University, and batch 2 (101 cases and 355 controls) was genotyped at the University of Michigan.1,155 cases and 1,117 controls were genotyped as part of the National Cancer Institute-sponsored Colorectal Transdisciplinary (CORECT) Study using a custom Affymetrix Axiom genome-wide platform measuring 1.2 million SNPs [48].3,768 cases and 3,028 controls were genotyped as part of the OncoArray Consortium using a custom Illumina OncoArray chip measuring 495 K SNPs (genome-wide backbone and known cancer susceptibility loci) [46].The final set of MECC participants including 173 cases and 140 controls were genotyped on the commercialized Infinium OncoArray-500 K Beadchip (Figure S1).
We imputed genotypes to the Haplotype Reference Consortium (HRC) panel [49] (39.2 million variants) using the University of Michigan Imputation Server [50], separately by genotyping platform.Imputed SNPs with allele count less than 20 and genotypes with quality scores (Rsq) less than 0.3 from any of the platforms were excluded from downstream analyses.Principal component analysis (PCA) was performed using PLINK 1.9 51 on directly genotyped SNPs shared across the four genotyping panels: Illumina Omni2.5, Affymetrix Axiom, Illumina Custom OncoArray, and Illumina Infinium OncoArray-500 K.After LD pruning (r 2 > 0.2), removing SNPs with minor allele frequency (MAF) < 0.01, and SNPs with PC1 and PC2 loading > 4.0, 55,852 autosomal SNPs were retained for PCA.Due to possible residual population substructure, the first 10 principal components for global ancestry were included in association analyses.

Quantification of tumor infiltrating lymphocytes
Hematoxylin and eosin-stained (H&E) tumor sections were blindly reviewed by a single gastrointestinal pathologist (J.K.G.) for 3,865 MECC participants.Tumor infiltrating lymphocytes (TILs) were evaluated as described previously [5].The mean count of TILs/high powered field (hpf ) for each tumor was calculated by taking the average number of TILs observed across five representative fields.Cases were separated into two groups (TILs/ hpf = 0 or TILs/hpf > 0) for analyses.

T cell receptor repertoire characterization
A subset of the 3,865 CRC cases that underwent pathology review had sufficient macrodissected tissue and extracted DNA available for immunosequencing.After QC, the analytic dataset included 2,750 cases.T cell receptor clonality and abundance data were derived from immunoSEQ assays run by Adaptive Biotechnologies (Seattle, WA).ImmunoSEQ uses a multiplex PCR system to amplify hypervariable complementarity determining region 3β (CDR3β) sequences of the TRB gene (T cell receptor beta locus; https://www.genenames.org/data/gene-symbol-report/#!/hgnc_id/HGNC:12155).Clonality was estimated by a modified Simpson diversity index which quantifies the clonality and diversity of the amino acid sequences of the T cell receptors.Abundance was estimated using the normalized number of TRB reads divided by the estimate of the total number of cells.We performed log transformation of clonality and abundance data due to the right-skewed distributions of the raw data.For clonality, we calculated a z-transformation of each sample average based on the distribution of those that underwent the same version of the assay.Two versions of the immunoSEQ assay were used, V2 (n = 1,165) and V4 (n = 1,585).Since the distributions of clonality were different between the two versions, all of the analyses of clonality were stratified by version and meta-analyzed.

Statistical analyses
Regression analysis using imputed genotype dosage was conducted using PLINK v1.9 [51] to investigate the associations between each SNP and each of the three T cell features (clonality, abundance, TILs).All models were adjusted for sex, age, genotyping center, and global ancestry.For clonality, analyses were performed separately by immunoSEQ assay version, and summary statistics were combined using fixed-effect inverse variance-weighted meta-analysis implemented in METAL [52].SNPs with a suggestive p value of less than 5 × 10 − 6 from the discovery phase were further assessed in the replication datasets and in a joint meta-analysis.

Replication phase Study populations and analyses
Three independent data sets were available for replication of our discovery phase findings: Colonomics ("CLX", N = 99; clonality, abundance, TILs/hpf ), the CRC Genetics Study ("CRCGEN", N = 162; clonality, abundance), and the Nurses' Health Study and the Health Professionals Follow-up Study (N = 505; TILs).The first replication dataset was a previously described set of 99 patients with microsatellite stable colon cancer diagnosed at stage II with immunoSEQ data.53, 54 The second replication data set (Colorectal Cancer Genetics Study; CRCGEN) included FFPE and frozen samples from 162 stage II microsatellite stable colon cancer patients with immunoSEQ data [54].
The third replication dataset was derived from two prospective cohort studies in the US: the Nurses' Health Study and the Health Professionals Follow-up Study [13,55].Tissue sections for all CRC cases were examined by a single pathologist (S.O.).TILs were dichotomized into absent and present [2,56].The second pathologist (J.N.G.) independently re-reviewed 398 randomly selected cases, and a good inter-observer correlation was observed, as previously described [2].

Imputation, QC, and statistical analyses for replication data sets
Genotypes from CLX and CRCGEN were imputed to the Haplotype Reference Consortium (HRC) panel [49] and underwent similar QC to the discovery phase.Immunosequencing data for clonality and abundance from both replication sets underwent the same QC measures and transformations as the MECC discovery set.Regression models using imputed genotypes were conducted to examine the associations with T cell features for CLX and CRCGEN studies, adjusted for sex, age, and global ancestry.
Genotype data from the Harvard cohorts were imputed to the 1000 Genomes Project panel and underwent QC steps as described previously [13].Logistic regression analyses were performed to investigate the association between each SNP and the presence of dichotomized TILs, after adjusting for sex, age, principal components, and tumor location.

Meta-analysis
A meta-analysis for the discovery and replication phases on the three immune-related outcomes using fixed-effect models with inverse variance weighting was implemented in METAL [52].Heterogeneity was evaluated using Cochran's Q test for heterogeneity and the measure I 2 .

Post-GWAS annotation
Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA GWAS, https://fuma.ctglab.nl/) was used for post-GWAS analyses for the SNPs with a p-value less than 5 × 10 − 6 in the discovery phase [57].eQTL mapping and chromatin interaction mapping were conducted using various data sources including Genotype-Tissue Expression (GTEx) and Hi-C data from FUMA GWAS.For SNPs with discovery p values < 5 × 10 − 6 for any of the three outcomes, additional eQTL analyses were performed as previously described using expression data on healthy colon mucosa samples (N = 485) from BarcUVa-Seq (https://barcuvaseq.org/) [58].We also examined the associations between the 19 SNPs with p values < 5 × 10 − 6 in joint meta-analyses and expression data from colorectal tumor tissues for genes located within 1 MB of the variants of interest using Colonomics (https://www.colonomics.org).[59]

Subset analyses within microsatellite stable (MSS) tumors
To examine whether the top variants identified from our study were associated with T cell features independent of MSI status, we performed the same analyses restricted to patients with MSS tumors.

Previously reported variants
We examined the associations between four germline variants associated with immune infiltration in tumors [8,9] as well as 155 variants from previous GWAS of CRC risk [11][12][13][14] and immune-related outcomes in our Discovery GWAS and joint meta-analyses.We also examined 1,587 variants associated with 33 immune traits with suggestive p < 10 − 6 from Sayaman et al. [9] in our discovery GWAS and joint meta-analyses.

Multiplex immunofluorescence analysis
To further investigate the associations between the most statistically significant variants identified from our GWAS and the density of specific tumor-associated T cell populations, multiplex immunofluorescence (mIF) was performed on tissue microarray (TMA) blocks containing tumor cores from 357 MECC participants [60].
TMA sections were immunostained using PerkinElmer OPAL™ 7-Color Automation Immunohistochemistry (IHC) kits (Waltham, MA) on the BOND RX autostainer (Leica Biosystems, Vista, CA).Panel markers included CD3, CD4, CD8, FOXP3, CD45RO (PTPRC), KRT (keratin), and DAPI to stain nuclei.Immunostained slides were imaged with the Vectra3® Automated Quantitative Pathology Imaging System.Multispectral TIFF images were exported for quantitative image analyses in HALO (Indica Labs, New Mexico), and cell densities for each fluorescent marker in the nucleus and in the stroma as well as proportion of cells positive for a given marker or marker combination were generated.Here, we specifically examined the densities (cells/mm 2 ) of cytotoxic T cells (CD3 + /CD8 + ), memory T cells (CD3 + /CD45RO + ), and regulatory T cell (CD3 + /CD4 + /FOXP3 + ) populations.We performed the inverse hyperbolic sine transformation on the T cell density data to account for skewed distributions.The density of regulatory T cells was dichotomized into two groups (0 or any) due to the heavy-zero distribution.
Imputed genotypes from the 19 top variants from our GWAS were used to examine the association between SNPs and the infiltration of cytotoxic, memory, and regulatory (0 vs. any) T cells as implemented in regression models adjusting for sex, age, genotyping center, and global ancestry.

Table 3
Two novel loci associated with TILs identified from Discovery GWAS and Meta-analyses with p < 5 × 10 − 6

Table 4
Descriptive statistics for discovery and replication datasets